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Abstract 

The  non-equilibrium  sediment  transport  (NEST)  modeling  approach  has  been  extended 
and  applied  in  this  chapter  to  simulate  non-cohesive  sediment  transport  induced  by 
rapidly-varying  transient  flows,  by  coastal  current  and  waves,  in  vegetated  water  bodies, 
and  by  overland  flow.  Even  though  different  flow  models  are  used  in  these  cases,  the 
sediment  transport  models  are  similar,  with  differences  in  sediment  entrainment, 
adaptation  length,  and  effective  diffusivity.  In  the  case  of  rapidly-varying  transient  flows, 
the  generalized  shallow  water  flow  equations  are  adopted  to  consider  interactions 
between  flow,  sediment  transport  and  bed  change.  In  the  coastal  context,  the  flow  model 
adopts  the  phase-averaged  shallow  water  flow  equations  with  wave-induced  radiation 
stresses  coupled  with  a  spectral  wave  transformation  model,  and  the  sediment  transport 
model  accounts  for  sediment  entrainment  and  mixing  (diffusion  and  dispersion)  by 
currents  and  waves.  In  the  case  of  vegetated  channels,  the  vegetation  drag  and  inertia 
forces  are  considered  in  the  momentum  equations  and  the  sediment  transport  capacity  is 
modified  due  to  vegetation  effect.  For  upland  soil  erosion,  the  overland  flow  is  simulated 
using  a  2-D  diffusion  wave  model  and  the  rill/interrill  erosion  due  to  raindrop  splash  and 
hydraulic  shear  is  considered  in  the  sediment  entrainment.  In  addition,  a  general  NEST 
model  framework  has  been  developed  for  simulating  transport  of  cohesive/non-cohesive 
sediment  mixtures,  taking  into  account  the  effects  of  cohesive  sediment  flocculation,  bed 
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consolidation  and  interactions  between  cohesive  and  non-cohesive  bed  materials. 
Selected  test  cases  demonstrate  that  the  extended  NEST  models  can  reasonably  reproduce 
the  sediment  transport  and  morphology  evolution  under  these  complex  flow  conditions. 

Keywords:  Cohesive  sediment;  non-cohesive  sediment;  rapidly-varying  transient  flows; 
upland  soil  erosion;  vegetation  effect;  waves 


Introduction 

To  mitigate  natural  and  man-caused  disasters  and  conserve  soil  and  water  resources, 
better  understanding  and  more  reliable  prediction  of  sediment  transport  by  various  surface 
water  flows  are  in  high  demand.  For  example,  rapidly-varying  transient  flows,  such  as 
dam/levee  break  flows,  storm  surge,  and  tsunami  waves,  cause  intensive  erosion  over  the 
surface  of  earthen  embankments  and  barriers  and  non-equilibrium  transport  of  non-uniform 
cohesive/non-cohesive  sediments  in  mixed  flow  regimes.  Soil  erosion  and  transport  by 
overland  flow  involve  raindrop  splash  and  rill/interrill/gully  erosion.  In  the  coastal 
environment,  sediment  transport  is  affected  by  both  currents  and  waves.  In  a  vegetated 
channel,  flow  is  retarded  and  diverted  by  vegetation  patches  and  induce  different  sediment 
transport  and  morphological  change  patterns,  whereas  in  a  coastal  marshland,  vegetation 
plays  an  important  role  in  stabilizing  the  bed  against  the  impacts  of  surges  and  waves.  In 
estuaries,  reservoirs,  lakes  and  rivers,  cohesive  and  non-cohesive  sediment  mixtures  widely 
exist  and  experience  complex  interactions  between  sediment  size  classes.  All  these  sediment 
transport  problems  are  quite  different  from  the  non-cohesive  sediment  transport  under 
common  channel  flow  conditions  described  in  the  previous  chapter  and  require  special 
modeling  capabilities.  Described  in  this  chapter  are  several  newly  developed  non-equilibrium 
sediment  transport  (NEST)  models  for  these  specific  applications. 


Sediment  Transport  under  Rapidly-Varying  Transient 

Flows 

The  morphodynamic  processes  under  actions  of  rapidly-varying  transient  flows,  such  as 
dam/levee  break  flows,  strong  storm  surge,  and  tsunami  waves,  involve  mixed 
supercritical/subcritical  flow  regimes,  hydraulic  jump,  strong  sediment  transport,  and  rapid 
bed  erosion.  Therefore,  simulation  of  flow,  sediment  transport,  and  morphology  evolution 
during  such  events  are  very  challenging.  To  accurately  simulate  these  processes,  the  flow 
model  is  required  to  have  shock-capturing  capabilities.  For  example,  the  1-D  and  depth- 
averaged  2-D  models  often  adopt  approximate  Riemann  solvers  and  TVD  (Total  Variation 
Diminishing)  schemes  (Toro,  2001),  while  the  vertical  2-D  and  3-D  models  use  the  VOF 
(volume -of-fluid)  method  and  smooth  particle  hydrodynamics  method  to  handle  the  rapidly 
varying  unsteady  flows.  Moreover,  the  model  should  consider  the  interactions  between  flow, 
sediment  transport,  and  bed  change.  There  have  been  several  numerical  models  in  the 
literature  to  simulate  earthen  embankment  break  flow  over  movable  beds  (e.g.,  Cao  et  al., 
2004;  Wang  and  Bowles,  2006;  Faeh,  2007;  Wu  and  Wang,  2007;  Roelvink  et  al.,  2009;  Wu 
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et  al.,  2012).  An  example  of  such  models  is  described  below,  which  improves  on  the  model  of 
Wuetal.  (2012). 

Model  Formulations 


When  the  sediment  concentration  is  high,  the  3-D  continuity  equation  of  water  and 
sediment  mixture  is  given  by  Equation  (1),  where  t  is  time,  xi  is  the  i'h  coordinate  (x,y,z 
for  i  =  1, 2, 3 ),  and  u  is  the  velocity  component  of  the  water  and  sediment  mixture  in  the  i,h 
direction.  In  Equation  (1),  p  is  the  combined  density  of  water  and  sediment  mixture  and  is 
given  by  p  =  pw(l  —  c)  +  psc ,  where  pw  and  p.  are  the  water  and  sediment  densities, 
respectively,  and  c  is  the  sediment  concentration  by  volume. 

dp  d(pni) 


=  0 


(1) 


dt  dxt 

By  using  the  Reynolds  time-averaging,  one  can  have  ui  =  uj  +  u\ ,  p  =  p  +  p  ,  c  =  c  +  c' , 
and  Equation  (1)  can  be  rewritten  as  Equation  (2).  In  Reynolds  time -averaging  the 
instantaneous  variables  (e.g.,  ui )  are  split  into  time-averaged  quantities  (e.g.,  uj )  and 
fluctuating  quantities  (e.g.,  u] ).  The  fluctuating  density,  p' ,  is  related  to  the  fluctuating 
sediment  concentration,  c' ,  by  p' =  (ps  -  pw)c' ,  so  that  p'li'  =  ( ps  —  pw)c'u[ .  The  turbulent 
sediment  flux  c'u'  is  modeled  by  using  the  sediment  turbulent  diffusion  concept,  yielding  the 
modeled  time-averaged  continuity  equation  as  given  by  Equation  (3),  where  es  is  the 
sediment  turbulent  diffusion  coefficient. 

dp  d(pu,)  d(p'u' ) 


dt  dxi 
dp  d(pUj) 


dX; 


=  0 


dt 


dx, 


=  (P,  “A)t- 
dx, 


dc A 
dx. 


(2) 


(3) 


'vi  V  '"'i  J 

Integrating  Equation  (3)  over  the  flow  depth  leads  to  the  depth-averaged  2-D  continuity 
equation  given  by  Equation  (4),  where  U  and  V  are  the  depth-averaged  flow  velocities  in  x 
and  y  directions  in  the  horizontal  plane,  h  is  the  flow  depth,  z,h  is  the  bed  elevation,  Ct  is  the 
depth-averaged  concentration  of  sediment  (total  load),  p  is  the  density  of  the  water  and 
sediment  mixture  in  the  water  column  determined  by  p  =  pw(\-C^)  +  psCt  (where  the 
overbars  are  omitted  for  simplicity).  The  density  of  the  water  and  sediment  mixture  in  the  bed 
surface  layer,  denoted  as  ph ,  is  determined  by  pb  =  pwp'm  +  ps  ( 1  —  p'm ) ,  where  p'm  is  the 

porosity  of  bed  material  in  the  surface  layer.  Since  Equation  (4)  governs  the  continuity  of  the 
water  and  sediment  mixture,  it  accounts  for  the  effect  of  sediment  turbulent  diffusion  by 
including  the  term  on  its  right-hand  side. 


d(ph)  d(pUh)  d(pVh)  dzb  , 

dt  dx  dy  "  dt  [  Ps  Pw’ 


— I  £,h 

dx 


dC, 

dx 
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Similarly,  depth-integrating  the  time-averaged  3-D  momentum  equations  of  the  water  and 
sediment  mixture,  using  the  hydrostatic  pressure  assumption,  and  ignoring  some  high-order 
turbulent  correlation  terms,  one  can  derive  the  depth-averaged  momentum  equations  of 
sediment-laden  shallow  water  flow  as  given  by  Equations  (5)  and  (6).  In  these  equations,  zs 
is  the  water  level,  n  is  the  Manning  roughness  coefficient,  g  is  the  gravitational  acceleration, 


2  +V2 


and  mb  = 


=  \ll  +  (dzb/dx)2  +  (dzb/dy)2 


/dt  is  the  eddy  viscosity,  U  =  y / 
considers  the  bed  slope  effect  but  can  be  lumped  in  the  Manning’s  n  coefficient. 

±(pUH)  +  l(pU^)  +  ±.{pUVh) 


dt 

df  .dU']  dr 


dx 


uth -  H - 

V  dx  J  dy 


i  dU 

v  dy  j 


dy 

dzs  1  , 2  dp  n2mbUU 

-psh^-ish'T,-pg^- 


which 


(5) 


f(^+£(pWq+|(prV) 


d  (  ,  dr 


dx 


H,h —  H - 

v  dx  J  dy 


dV 


v  dy , 


\  ,  dzs  1,2  dp  n2mbUV 

f  "‘V?  ^y-pg-fF- 


(6) 


The  sediment  transport  and  bed  change  are  determined  by  Equations  (7)  and  (8),  where 
qt,  is  the  capacity  or  equilibrium  transport  rate  of  total-load  (bed-material  load)  sediment, 

and  Lt  is  the  adaptation  length  of  sediment  (see  the  previous  chapter).  Note  that  the 
coefficient  mb  considers  the  lateral  erosion  along  the  side  slope.  It  might  be  lumped  in  the 
sediment  transport  capacity  qtt .  For  many  rapidly-varying  transient  flows,  the  turbulent  stress 
and  diffusion  terms  in  Equations  (4)— (7)  are  usually  ignored  because  they  are  much  smaller 
than  the  convection  terms. 


8{hC, )  d(hUCt )  8(hVCt ) 


dt 


dx 


dy 


d(  dC, 
s  h — - 
dx 


dx 


+  d_ 
dy 


'sh8^ 

vs  dy  j 


y(ukC,-mM±)  (7) 


■ 


(8) 


The  sediment  transport  capacity  qlt  is  determined  using  the  formula  of  Wu  et  al.  (2000), 
with  the  modification  suggested  by  Wu  (2007)  to  consider  the  effect  of  gravity  on  sediment 
transport  over  steep  slope  as  explained  in  the  previous  chapter.  In  addition,  to  consider  the 
effect  of  sediment  concentration,  the  settling  velocity  cos  used  in  Wu  et  al.  (2000)  formula  is 


determined  by  Richardson  and  Zaki’s  (1954)  method  as  cos  =<ys0(l  —  C,)"',  with  m  «  4  and 
a>s0  being  the  settling  velocity  of  single  sediment  particles  in  quiescent,  distilled  water. 

Considering  nearly  vertical  sidewalls  or  banks  observed  in  laboratory  and  field 
experiments  even  for  non-cohesive  sediment  materials,  a  specially  designed  avalanching 
algorithm  is  adopted  here.  The  non-cohesive  slope  avalanching  algorithm  presented  by  Wu 
(2007)  is  modified  by  specifying  two  different  repose  angles  for  the  submerged  materials  and 
the  emergent  materials  above  the  water  surface. 

By  using  Equations  (7)  and  (8),  Equations  (4) — (6)  can  be  reformulated  as  Equations  (9)— 
(11),  respectively.  The  reformulated  continuity  equation,  Equation  (9),  is  much  simpler  than 
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Equation  (4).  Moreover,  the  flow  density  is  eliminated  from  the  left-hand  sides  of  Equations 
(9)— (1 1),  so  that  the  traditional  shock-capturing  schemes  can  be  more  conveniently  adopted  to 
solve  these  equations. 

dll  |  d(Uh)  [  d(Uh)  _  1  UhCt  —  mbqtt 

dt  dx  dy  1  -p'm  Lt 

d(Uh )  d(u2h)  d(UVh) 
dt  dx  dy 

1  df  .  dU  \  df  ,dU)  dzs  1  , 2  1  dp  n2mbUU 

~~p\dx\!l,h ~dx)  +  dy[Mth~dyj\~8h  ~dx  ~2gh  ~p~dx~  8  hm 

+u^\-^{s  Afg 

p  [  dx\s  &  )  5y^  *  dy  J  ^  1  -p'mJ  L, 

d(Vh)  d(UVh)  d(v2h) 
dt  dx  dy 

if  df  ,5V df  ,5V ^l]  dz.  1  l  dp  n2m.UV 

p  dx{  dx  J  dyy  dy  J  dy  2  p  dy  h13 


■  Ps  Pw 


- £.h 


- £,h 


UhC,  -  mbq, , 

L, 


Wu  et  al.  (2012)  solved  Equations  (9)— (11)  without  the  turbulent  diffusion  terms  using  an 
explicit  finite  volume  method  on  a  rectangular  grid.  The  HLL  approximate  Riemann  solver 
(Harten  et  al.,  1983)  is  used  to  determine  the  convective  fluxes  at  cell  faces  in  both  x  and  y 
directions.  A  second-order  accurate  approximation  in  space  is  obtained  by  reconstructing  the 
left  and  right  states  or  fluxes  using  the  Monotonic  Upstream  Scheme  for  Conservation  Laws 
(MUSCL)  (van  Leer,  1979).  In  a  later  study,  Wu  and  Marsooli  (2012)  adopted  a  hybrid 
approach  in  which  the  streamwise  convective  flux  is  treated  using  the  HLL  solver  and 
MUSCL  reconstruction,  but  the  lateral  convective  flux  is  determined  using  the  HLPA  (Hybrid 
Linear/Parabolic  Approximation)  scheme  (Zhu,  1991),  which  has  approximately  second-order 
accuracy  in  space.  This  hybrid  approach  has  less  numerical  diffusion  than  the  former 
approach  applying  the  HLL  solver  in  both  x  and  y  directions. 

The  water-surface  gradient  terms  in  Equations  (10)  and  (11)  are  treated  using  the 
modified  center  difference  scheme  suggested  by  Ying  and  Wang  (2008).  If  the  turbulent 
diffusion  terms  are  included,  they  can  be  discretized  using  the  traditional  center  difference 
scheme.  The  convective  flux  of  sediment  at  cell  faces  is  determined  using  the  HLPA  scheme. 
The  bed  change  equation,  Equation  (8),  is  discretized  using  a  point-wise  forward  difference 
scheme.  It  is  well  known  that  oscillations  in  flow  velocity  near  hydraulic  jumps  may  be 
generated  by  the  HLL  approximate  Riemann  solver.  In  order  to  eliminate  the  effects  of  these 
numerical  oscillations  on  sediment  transport  and  bed  changes,  the  velocity  used  in  the 
calculation  of  sediment  transport  capacity  is  set  as  the  convective  flux  averaged  from  the  two 
cell  faces  in  x  or  y  direction  divided  by  the  flow  depth  (Wu  et  al.,  2012). 
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The  model  was  tested  using  a  benchmark  experiment  carried  out  at  Universite  catholique 
de  Louvain  (UCL),  Belgium  to  investigate  the  two-dimensional  morphological  evolution  of  a 
movable  bed  under  the  action  of  partial  dam-break  wave  (IAHR  Working  Group,  2012).  The 
flume  was  3.6  m  wide  and  about  36  m  long,  as  shown  in  Figure  1.  The  partial  dam  break  was 
represented  by  rapidly  lifting  the  1-m  wide  gate  between  two  impervious  blocks.  The  fixed 
bed  of  the  flume  was  covered  with  an  85  mm  thick  sand  layer  extended  from  1  m  upstream  of 
the  dam  to  9  m  downstream  of  the  gate.  The  sand  was  fully  saturated  before  the  experiment 
and  had  a  median  diameter  of  1.61  mm,  a  specific  gravity  of  2.63,  and  an  initial  bed  porosity 
of  0.42.  The  initial  water  level  in  the  reservoir  was  0.47  m  above  the  fixed  flume  bed,  while 
the  downstream  was  dry  bed  in  this  test  case.  A  mesh  of  600x140  cells  in  longitudinal  and 
lateral  directions  was  used  to  cover  the  entire  flume.  The  longitudinal  and  lateral  grid  spacing 
was  0.025  m  in  the  movable  bed  reach  near  the  breach  width  and  increased  gradually  in  the 
downstream  and  upstream  directions  and  towards  both  sidewalls.  A  constant  time  step  of 
0.0005  s  was  used  in  this  case.  The  Manning’s  n  in  the  sand  bed  was  suggested  as  0.0165  by 
the  benchmark  test  (IAHR  Working  Group,  2012).  For  the  portion  of  fixed  bed,  the 
Manning’s  n  was  set  as  0.01.  The  suspended -load  adaptation  coefficient  a  was  set  as  4.0 
and  the  bed-load  adaptation  length  Lb  was  0.025  m. 


Figure  1 .  Configurations  of  the  UCL  partial  dam  break  experiment  (dimension  in  meters) 


Figures  2(a)  and  (b)  compare  the  calculated  and  measured  bed  topographies  at  the  end  of 
the  experiment  ( t  =  20  s).  The  measured  bed  profiles  along  three  longitudinal  sections 
( y  =  0.2 ,  0.7  and  1.45  m)  at  f  =  20  s  are  compared  with  the  calculated  results  in  Figure  3. 
Note  that  Figure  3  includes  bed  profiles  measured  in  two  experimental  runs  using  the  same 
setup,  and  difference  between  the  two  measurements  shows  the  difficulty  of  repeating  such 
experiments.  Both  calculation  and  measurement  showed  that  significant  erosion  led  to  a  large 
scour  hole  downstream  of  the  dam.  Some  of  the  eroded  sediment  deposited  along  the 
sidewalls  where  hydraulic  jumps  appeared  and  downstream  of  the  scour  hole.  The  deposition 
pattern  was  controlled  by  the  cross  waves  generated  by  the  partial  dam  break,  which  spread  to 
the  sidewalls  of  the  flume  and  were  then  reflected  from  one  side  to  the  other  side.  The 
calculated  scour  hole  was  deeper  than  the  measured  one. 
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Figure  2.  Bed  topographies  at  r=20  s:  (a)  calculated  and  (b)  measured 
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Figure  3.  Measured  and  calculated  final  bed  profiles  along  three  longitudinal  sections 
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Figure  4  shows  the  measured  and  simulated  time  series  of  water  surface  elevations 
(above  the  fixed  flume  bottom)  at  8  gage  points  marked  in  Figure  1.  The  variations  of  water 
surface  at  Gages  1  and  4  have  the  same  trend,  due  to  almost  symmetric  flow,  sediment 
transport,  and  bed  change  patterns  about  the  channel  centerline  observed  by  the  model 
simulation.  Similar  trends  can  also  be  observed  between  Gages  2  and  3,  between  Gages  5  and 
8,  and  between  Gages  6  and  7.  The  calculated  water  levels  at  Gages  1  and  4  are  lower  than 
the  measured  ones.  This  may  be  due  to  strong  3-D  features  of  the  flow  near  the  two  corners  of 
sudden  expansion  of  the  flume  width,  as  well  as  the  over  predicted  erosion  there. 


Figure  4.  Measured  and  calculated  water  surface  variations  with  time  at  8  gage  points 


The  model  was  also  tested  using  the  IMPACT  laboratory  experiments  carried  out  by  HR 
Wallingford,  UK  (Morris  et  al.,  2009).  The  flume  was  50  m  long  and  10  m  wide.  A  model 
dam  of  0.5  m  in  height  was  located  at  about  36  m  downstream  of  the  flume  entrance.  It  had  a 
crest  width  of  0.2  m  and  up/downstream  slopes  of  1V:1.7H.  A  wide  grading  and  a  narrow 
grading  sediment  mixture  were  tested  in  the  experiment.  The  two  sediment  mixture 
gradations  had  the  same  median  diameter  of  0.25  mm.  A  uniform  mesh  consisting  of 
1000x200  nodes  in  longitudinal  and  lateral  directions  was  used  in  the  simulation.  The  dam 
material  porosity  was  estimated  as  0.36,  and  the  Manning’s  n  was  estimated  as  0.018.  The 
initial  pilot  channel  was  0.02  m  deep  and  0.152  m  wide.  The  repose  angle  of  the  bed  material 
was  set  to  be  31°  in  the  submerged  condition  and  about  80°  in  the  emergent  condition.  The 
inflow  discharge  of  0.07  m 7s  was  specified  at  the  flume  entrance.  Calibration  showed  that  a 
constant  Lt  of  0.07  m  gave  reasonable  results  in  this  case. 

Figure  5  shows  the  calculated  bed  topography  and  flow  pattern  at  the  peak  breach  flow. 
At  the  downstream  of  the  breach,  the  flow  spread  laterally  and  small  bed  forms  or  undulations 
appeared  on  both  sides  of  the  flow.  The  simulated  breach  exhibits  a  “bell  mouth”  weir  flow 
and  topography,  as  observed  in  the  experiment.  Figure  6  compares  the  measured  and 
calculated  breach  flow  discharges  variation  with  time,  and  Figure  7  compares  the  measured 
and  calculated  breach  widths  and  headwater  levels.  The  two  figures  show  that  the  model 
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predicted  a  steeper  rising  limb  of  breach  flow  hydrograph,  a  slower  tailing  of  reservoir  water 
level  and  a  faster  widening  process  than  the  experiment.  The  calculation  results  agree 
generally  well  with  the  measured  data.  In  particular,  the  peak  breach  flow  discharge  and  the 
final  breach  width  are  well  predicted  by  the  model. 


Figure  5.  Simulated  flow  pattern  and  bed  topography  at  the  peak  breach  flow 


Figure  6.  Measured  and  simulated  breach  flow  discharges 


Figure  7.  Measured  and  simulated  breach  widths  and  headwater  levels 
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Sediment  Transport  by  Currents  and  Waves 

In  the  context  of  coastal  environment,  because  the  multiple  interacting  forces  (waves, 
wind,  tide,  river  flow,  density  current,  etc.)  exhibit  a  wide  range  of  spatial  and  temporal 
scales,  the  complex  costal  physical  processes  are  quantitatively  not  well  understood. 
Traditionally,  most  coastal  sediment  transport  models  adopt  the  equilibrium  transport 
concept.  The  following  section  presents  an  application  of  the  NEST  model  in  the  Coastal 
Modeling  System  (CMS)  as  illustrated  by  Sanchez  and  Wu  (201 1). 

Model  Formulations 

The  CMS,  developed  under  the  U.  S.  Army  Corps  of  Engineers’  Coastal  Inlets  Research 
Program  (CIRP),  is  designed  for  practical  applications  in  navigation  channel  performance 
evaluation  and  sediment  management  near  coastal  inlets  and  adjacent  beaches.  It  is  a  phase- 
averaged  model,  in  which  the  circulation  model  computes  the  unsteady  water  level  and 
current  velocity  fields  by  solving  the  phase-averaged,  depth-averaged  2-D  shallow  water  flow 
equations  on  a  non-uniform  or  quadtree  Cartesian  grid  with  an  explicit  and  implicit  finite 
volume  scheme  (Buttolph  et  al.,  2006;  Wu  et  al.,  2011).  The  model  can  simulate  tide,  wind 
and  wave  driven  currents,  and  includes  the  Coriolis  force,  wind  stress,  bottom  friction,  and 
wave  radiation  stresses.  The  phase-averaged  2-D  shallow  water  flow  equations  can  be  written 
as  Equations  (11)— (13)  of  the  previous  chapter.  In  these  equations  Z>  ( i,  j  =  x,y  )  represents 

the  wave-driven  radiation  stresses  and  is  given  by  Equation  (12),  where  E  is  the  wave  energy 
density,  /  is  the  wave  frequency,  6  is  the  wave  direction,  wt  is  the  wave  unit  vector 

(cos 6 ,  sin 1 9) ,  ng  =  0.5(1  +  2/c/?/sinh  2kh ) ,  and  k  is  the  wave  number.  The  bed  shear  stresses 

experienced  by  the  currents  in  the  case  where  both  waves  and  currents  coexist  are  given  by 
Equations  (13)  and  (14),  where  Uwm  is  the  maximum  orbital  bottom  velocity  of  waves  and 

cf  =gn2/hm  . 

Dij  =  JJ E(f  ,0)\jigwiwj  +  Sv  (ng  -0.5 JjdfdO  (12) 

Tbx=cfp{u2+V2+0.5Uin)l'2U  (13) 

Tby=cfp(u2 +v2  +  0.5uln)ll2V  (14) 

The  wave  characteristics  are  determined  using  a  spectral  wave  transformation  model  that 
solves  the  steady-state  wave-action  balance  equation  on  a  non-uniform  Cartesian  grid  with  a 
finite  difference  scheme.  It  considers  wind  wave  generation  and  growth,  diffraction, 
reflection,  dissipation  due  to  bottom  friction,  white  capping  and  breaking,  wave-wave  and 
wave -current  interactions,  wave  runup,  wave  setup,  and  wave  transmission  through  structures 
(Lin  et  al.,  2008).  The  sediment  transport  model  calculates  the  depth-averaged  2-D  non¬ 
equilibrium  transport  of  total  load  and  the  resulting  bed  change  using  Equations  (15)  and  (16), 
where  rs  is  the  ratio  of  suspended  load  to  total  load,  U c  is  the  current  speed,  D  is  the 
empirical  coefficient  for  the  effect  of  bed  slope  on  bed  load,  and  the  subscript  j  denotes  the 
index  of  horizontal  coordinates  x  and  y  . 
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Note  that  the  effect  of  bed  slope  on  bed-load  transport  is  taken  into  account  by  the  last 
term  in  Equation  (16),  whereas  this  is  accounted  for  through  the  transport  direction  cosines 
abx  and  ah  in  Equation  (24)  in  the  previous  chapter.  The  approach  of  Equation  (16)  is  found 


to  be  more  stable,  but  with  a  loosely  defined  coefficient  Ds .  Watanabe  (1985)  used  D  =  10 
based  on  a  free-body  diagram  of  the  grain  in  a  sloping  bed.  Struiksma  et  al.  (1985)  used 
=  4.  Later  studies,  such  as  Larson  et  al.  (2003)  and  Karambas  (2003),  reported  good 


results  with  Ds  =  2.  Sanchez  and  Wu  (2011)  used  a  value  of  about  1  to  avoid  over-smoothing 
the  morphology  change.  In  practice,  Ds  may  be  a  function  of  the  flow  and  sediment 
characteristics  and  vary  from  site  to  site. 

The  effective  sediment  diffusivity  es  is  related  to  the  eddy  viscosity  vt  by  s%  =  vi  / cr  ,  in 


which  cr  is  the  Schmidt  number.  The  eddy  viscosity  includes  the  effects  of  currents,  waves, 
and  wave  breaking  as  given  by  Equation  (17),  where  vtc  is  the  eddy  viscosity  due  to  currents 
and  can  be  determined  using  several  turbulence  models,  such  as  the  depth-averaged  parabolic 
model  and  the  modified  mixing-length  model,  A  is  an  empirical  coefficient  representing  the 
lateral  wave  mixing  strength,  H  is  the  wave  height,  Db  is  the  energy  dissipation  rate  due  to 


wave  breaking,  and  kb  is  an  empirical  coefficient.  The  second  and  third  terms  on  the  right- 
hand  side  of  Equation  (17)  account  for  the  effects  of  waves  and  wave  breaking,  respectively. 


v,  =v,+AU . H+t 


A,  , 

V  P 


(17) 


The  sediment  transport  capacity  Ct,k  under  currents  and  waves  is  significantly  affected 
by  interactions  between  currents  and  waves.  The  available  formulas  for  determining  Cf>k 
include  Watanabe  (1987),  Soulsby  (1997),  Lund-CIRP  (Camenen  and  Larson  2007),  and  van 
Rijn  (2007a,  b).  Wu  and  Lin  (2012)  extended  the  non-uniform  sediment  transport  capacity 
formula  of  Wu  et  al.  (2000)  to  the  cases  where  both  waves  and  currents  exist. 


Model  Test  Example 

The  model  was  tested  against  many  laboratory  and  field  case  studies  (Sanchez  and  Wu, 
2011;  Sanchez  et  al.,  2011).  One  of  the  field  cases  was  the  application  in  the  Shinnecock  Inlet 
on  Long  Island,  NY.  It  is  a  mixed-energy,  wave  dominated  inlet,  with  mainly  semidiurnal  tide 
with  an  average  spring  tidal  range  of  1.1  m.  The  Shinnecock  Bay  has  water  depths  typically 
less  than  2  m  and  a  tidal  prism  of  3.29x0  m  .  Typical  wave  conditions  are  from  the  southeast 
with  wave  heights  of  1  m  and  periods  of  7  s,  while  northeast  storms  and  hurricanes  can 
produce  wave  heights  in  excess  of  4  m  with  periods  of  12—14  s  (Buonaiuto  and  Militello, 
2003).  The  simulation  covered  the  time  period  between  August  13,  1997  and  May  28,  1998. 
Figure  8  shows  the  computational  domain  with  initial  bathymetry  of  1997.  The  grid  for  flow 
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computation  had  49,780  cells  and  a  variable  grid  resolution  between  15  and  100  m.  The 
sediment  size  ranged  from  0.125  mm  offshore  and  inside  the  bay  to  2  mm  in  the  inlet,  but  the 
predominant  size  was  0.5  mm,  which  was  used  in  the  simulation.  The  sediment  transport 
capacity  was  calculated  with  the  Watanabe  (1987)  equation  adjusted  by  multiplying  a  factor 
of  0.25.  The  sediment  adaptation  length  was  given  as  20  m,  and  the  bed  slope  coefficient  Ds 


was  1.0. 

The  calculated  longshore  sediment  transport  rate  was  133,000  yd /yr,  which  was  within 
the  range  of  110,000-190,000  yd3/yr  reported  by  Rosati  et  al.  (1999).  Figure  9  shows  the 
measured  and  computed  morphology  changes.  Qualitatively,  the  results  showed  several 
features  in  common  such  as  deposition  on  the  peripheral  of  the  ebb  shoal,  erosion  at  both  the 
east  and  west  bypass  bars,  erosion  of  the  eastern  portion  of  the  south  bypass  bar,  and 
accretion  in  the  deposition  basin.  Quantitatively,  the  morphodynamic  model  correctly 
predicted  either  an  erosional  or  depositional  trend  at  approximately  66%  of  the  computational 
cells  (Sanchez  and  Wu,  201 1). 


Figure  8.  Computational  domain  with  initial  bathymetry  for  1997  (square  cells  represent 

inactive  land  cells) 


Figure  9.  (left)  Measured  and  (right)  simulated  morphology  changes  between  August  1997 
and  May  1998  (Sanchez  and  Wu,  201 1) 
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Sediment  Transport  in  Vegetated  Channels 


Governing  Equations 


Flow  through  vegetation  usually  is  unsteady  and  turbulent,  but  the  mostly  considered 
flow  properties  in  practice  are  the  time-  and  space-averaged  behaviors  rather  than  the  detailed 
instantaneous  flow  features  around  individual  vegetation.  By  double-averaging  the  Navier- 
Stokes  equations  in  time  and  space,  one  can  derive  the  governing  equations  of  free-surface 
flow  through  vegetation  (Nikora  et  al.,  2007;  Lopez  and  Garcia,  2001).  Based  on  the 
hydrostatic  pressure  assumption,  further  depth-averaging  (Wu  and  Wang,  2004a)  the  double- 
averaged  3-D  Navier-Stokes  equations  leads  to  the  shallow  water  equations  with  vegetation 
effects  as  given  by  Equations  (18)— (20),  where  cv  is  the  depth-averaged  volumetric 
concentration  of  vegetation,  and  Fdx  and  Fdy  are  the  x  -  and  y  -components  of  the  drag  force 
and  inertia  force  on  vegetation  exerted  by  the  flow.  These  forces,  as  defined  by  Morison  et  al. 
(1950),  are  given  by  Equation  (21),  where  Cd  is  the  drag  coefficient  of  vegetation,  CM  is  the 
inertia  coefficient  of  about  2,  p  is  the  water  density,  Nv  is  the  vegetation  density  defined  as 


number  of  vegetation  elements  per  unit  horizontal  (bed)  area,  A  is  the  projected  area  defined 
as  the  frontal  area  of  a  vegetation  element  projected  to  the  plane  normal  to  the  stream-wise 
flow  direction,  Vv  is  the  volume  of  a  vegetation  element,  Uv  is  the  vector  of  flow  velocity 
acting  on  the  vegetation  element,  and  | Ur  |  (or  simply  written  as  JJr )  is  the  magnitude  of  JJV . 


For  emergent  vegetation,  Uv  is  the  depth-averaged  flow  velocity  U  ,  while  for  submerged 
vegetation,  Uv  should  be  the  velocity  averaged  only  over  the  vegetation  layer,  as  shown  in 
Figure  10. 


d[(l-c>]  d[(l -cv)hU]  5[(1  -cv)hV] 
dt  dx  dy 


(18) 


a[(i -cv)hu]  s[(i -cr)hu2]  a[(i -cr)huv] 

- S, - + - T, - + - dy - +  s(W'),‘ 

1  a[(l-c,)/.7-,]  i  8[(1-C,.)/T,] 

p  dx  p  dy  p  p  dx 

d[{l-cv)hV]  |  d[{l-cv)hUV]  |  a[(l-c„)/zV2]  |  c 

dt  +  dx  +  dy  +8^  Cv’\ 
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The  depth-averaged  2-D  transport  equation  of  suspended  load  in  vegetated  channels  is 
given  by  Equation  (22),  where  Ck  is  the  depth-averaged  concentration  of  the  kth  size  class  of 


suspended  load.  The  vegetation  density  term  1  — cv  also  appears  in  the  double-averaged  3-D 
Navier-Stokes  equations  and  suspended-load  transport  equation  (Wu,  2007).  The  transport 
equation  of  bed  load  in  vegetated  channels  is  given  by  Equation  (23),  where  cvb  is  the 


volumetric  concentration  of  vegetation  in  the  bed-load  layer. 

d[(l -cv)hCk/j3sk]  |  d[(l -cv)hUCt]  |  d[(l-cv)hVCk] 
dt  dx  dy 

=  4-  (]-c,)6'Jld§^-  +  -c,) 

dx\_  ox  \  dy\_  dy 


S[{l~cvb)qbk/ubk]  d\abl{\-cvb)qbk~\  d[aby(\-cvb)qbk']  i 

- 3, - + - T, - + — - =  Z(1 ~c“ )(,J-  "  "*  > (23) 


The  bed  change  can  be  determined  by  Equation  (24),  where  cvbei  is  the  volumetric 
concentration  of  vegetation  at  the  bed  surface. 

V 

dtJk 

In  the  case  of  low  vegetation  density,  Equations  (18)— (20),  (22)  and  (23)  reduce  to  the  2- 
D  flow  and  sediment  transport  equations  expressed  in  Equations  (11)— (13),  (15)  and  (24)  in 
the  previous  chapter,  with  the  drag  force  and  inertia  force  added  in  the  momentum  equations. 

The  drag  coefficient  Cd  is  the  key  parameter  in  the  vegetation  effect  model.  It  is  related 
to  vegetation  properties  and  flow  conditions.  The  turbulent  stresses  can  be  determined  using 
the  modified  mixing-length  model  (Wu  and  Wang,  2004a)  and  k-s  turbulence  model 
(Lopez  and  Garcia,  2001).  The  turbulence  generated  by  vegetation  needs  to  be  considered  in 
the  k  and  s  equations.  For  details  refer  to  Wu  (2007). 


(1  -  cv )  (Ctk  -  C, )  +  -  (1  -  cvb )  (qh,k  -  qbk )  (24) 


0  -fC)0  -Gfcd) 


Vegetation  Control  volume 


(a)  Emergent  (b)  Submerged 

Figure  10.  Vegetation  in  open  channels  (side  view) 


Sediment  Transport  Capacity  Formulas 

Sediment  transport  capacity  in  vegetated  channels  is  affected  by  vegetation.  Okabe  et  al. 
(1997),  Jordanova  and  James  (2003),  and  Wu  et  al.  (2005)  found  that  the  bed-load  transport  is 
mainly  related  to  the  bed  shear  rather  than  the  drag  force  exerted  on  vegetation  elements.  If 
the  effective  bed  shear  is  used,  some  existing  empirical  formulas  developed  for  the  bed-load 
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transport  in  non-vegetated  channels  can  be  extended  to  vegetated  channels.  Jordanova  and 
James  (2003)  experimentally  investigated  the  bed-load  transport  in  a  flume  covered  with 
uniformly  distributed  vegetation.  They  used  the  method  of  Li  and  Shen  (1973)  to  determine 
the  effective  bed  shear,  and  proposed  a  simple  empirical  formula.  Okabe  et  al.  (1997)  used  the 
k—e  turbulence  model  to  compute  the  effective  bed  shear  stress,  and  found  that  the  Ashida- 
Michiue  (1972)  formula  can  be  used  to  determine  the  bed-load  transport  rate  in  vegetated 
channels.  Wu  et  al.  (2005)  and  Wu  and  He  (2009)  found  that  the  formula  of  Wu  et  al.  (2000) 
is  applicable  for  the  bed-load  transport  in  vegetated  channels  if  the  approach  of  Barfield  et  al. 
(1979)  is  used  to  determine  the  effective  bed  shear  stress  as  given  by  Equation  (25),  where  S 
is  the  bed  slope  or  energy  slope,  y  is  the  specific  weight  of  water,  and  Rs  =  h  ln j ( 2h  +  / , ) ,  in 
which  h  is  the  flow  depth,  and  ln  is  the  lateral  spacing  between  vegetation  elements. 

rh=yRS  (25) 

The  suspended-load  transport  in  vegetated  channels  has  not  been  extensively 
investigated.  As  vegetation  may  reduce  the  mean  flow  velocity  significantly  but  intensify  the 
turbulence  in  a  vegetated  zone,  the  effect  of  vegetation  on  suspended-load  transport  is  more 
complex.  More  experimental  and  theoretical  studies  are  needed  to  quantify  this  effect.  As  an 
approximation,  one  may  apply  some  existing  suspended  load  formulas  established  under  non- 
vegetated  conditions  to  vegetated  channels.  Certainly,  this  application  should  be  carried  out 
with  caution. 

Model  Test  Examples 

The  depth-averaged  2-D  model  of  rapidly-varying  transient  flows  described  in  the 
previous  section  has  been  extended  to  flow  in  vegetated  waters  by  adding  the  drag  and  inertia 
forces  in  the  momentum  equations,  Equations  (10)  and  (11).  The  model  has  been  tested  using 
Bennett  and  Alonso’s  (2003)  experiments  on  morphodynamics  near  alternate  vegetation 
zones  in  an  open  channel.  Two  experiment  runs  were  conducted  in  a  laboratory  flume  that 
had  a  preformed  trapezoidal  cross-section  with  side  slopes  of  34°,  and  two  alternate 
vegetation  zones  separated  by  a  distance  of  1.5  m,  as  shown  in  Figure  11.  The  vegetation 
zones  were  0.49  m  long  and  0.245  m  wide  rectangles  in  experimental  run  1  and  semi-circles 
with  radii  of  0.245  m  in  run  2.  In  both  experimental  runs,  the  vegetation  elements  were 
emergent,  rigid,  and  in  a  staggered  pattern,  and  had  a  stem  diameter  of  4.8  mm  and  a  density 
of  2500  units/m".  The  sediment  size  was  0.8  mm.  The  computational  mesh  is  uniform,  with 
Ax  =  0.041  m  and  Ay  =  0.016  m,  as  shown  in  Figure  11.  The  Manning’s  roughness 
coefficient  is  adjusted  to  0.028,  which  yields  a  uniform  flow  in  the  absence  of  vegetation.  The 
drag  coefficient  Cd  is  set  as  2.0.  The  sediment  adaptation  length  is  determined  using 

Equation  (32)  in  the  previous  chapter,  with  Lh  =  0.025  m  and  a  =  4  . 

In  this  test  case,  the  effects  of  sediment  concentration  and  bed  changes  on  the  flow  can  be 
ignored  because  the  flow  is  subcritical  and  the  sediment  transport  is  weak.  The  bank  slope  is 
high,  so  that  the  effect  of  gravity  on  bed-load  transport  on  the  slope  needs  to  be  considered. 
This  is  done  using  Equation  (16),  with  Ds=  10 .  As  the  sediment  particles  were  non-cohesive, 
the  model  simulates  the  bank  erosion  and  retreatment  by  applying  repose  angle  to  the  bank 
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slope.  The  wetted  and  dry  repose  angles  are  set  as  34°  and  50°,  respectively,  as  explained  in 
the  previous  chapter. 

Figures  12  and  13  show  the  comparison  of  the  measured  and  predicted  changes  in  bed 
and  bank  topography  near  the  second  vegetation  zone  in  an  elapsed  time  of  110  minutes. 
Deposition  occurs  in  and  near  the  vegetation  zone,  whereas  erosion  occurs  in  the  area  away 
from  the  vegetation  zones,  resulting  in  a  meandering  channel  planform.  The  channel  erosion 
and  deposition  magnitudes  and  distributions  are  predicted  reasonably  well.  The  erosion  and 
deposition  magnitudes  near  the  rectangular  vegetation  zone  are  larger  than  those  near  the 
semi-circular  vegetation  zone. 


Figure  11.  Mesh  and  sketch  of  the  experiments  of  Bennett  and  Alonso  (2003) 


00 


Figure  12.  Changes  in  bed  topography  due  to  rectangular  vegetation  zones:  (a)  Measured  by 
Bennett  and  Alonso  (2003);  (b)  Simulated.  (Dashed  lines  show  the  location  of  the  channel  top 
width,  and  rectangles  show  the  extent  of  the  vegetation  zone) 
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(a) 


Figure  13.  Changes  in  bed  topography  due  to  semi-circular  vegetation  zones:  (a)  Measured  by 
Bennett  and  Alonso  (2003);  (b)  Simulated.  (Dashed  lines  show  the  location  of  the  channel 
top  width,  and  semi-circles  show  the  extent  of  the  vegetation  zone) 


Soil  Erosion  and  Transport  by  Upland  Flow 

To  investigate  the  complex  hydrological,  morphodynamic,  and  environmental  processes 
in  watersheds,  a  physically-based  integrated  2-D  surface  and  3-D  subsurface  model  for  flow, 
soil  erosion  and  transport,  and  contaminant  transport  was  developed  by  He  et  al.  (2009).  The 
model  simulates  the  rainfall-induced  surface  flow  by  solving  the  depth-averaged  2-D 
diffusion  wave  equation  and  the  variably  saturated  subsurface  flow  by  solving  the  3-D  mixed- 
form  Richards  equation.  The  surface  and  subsurface  flow  equations  are  coupled  using  the 
continuity  conditions  of  pressure  and  exchange  flux  at  the  ground  surface.  The  NEST  model 
is  applied  to  compute  the  non-uniform  soil  erosion  and  transport  in  upland  fields,  considering 
detachments  by  rain  splash  and  hydraulic  erosion  driven  by  surface  flow.  The  integrated  2-D 
surface  and  3-D  subsurface  contaminant  transport  model  takes  into  account  the  contaminant 
changes  due  to  sediment  sorption  and  desorption,  as  well  as  exchanges  between  surface  and 
subsurface  domains  due  to  infiltration,  diffusion,  and  bed  change.  The  2-D  soil  erosion  and 
transport  model  component  is  briefly  presented  here. 


Model  Formulations 
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When  inertia  is  not  important,  the  surface  flow  can  be  described  using  the  depth-averaged 
2-D  diffusion  wave  equation  as  suggested  by  Di  Giammarco  et  al.  (1996)  and  given  by 
Equation  (26),  where  Hs  is  the  water  surface  elevation  (m),  which  is  the  sum  of  surface  water 
depth  h  (m)  and  bed  elevation  zb  (m),  qr  is  the  rainfall  rate  (m/s),  qe  is  the  water  exchange 
rate  with  subsurface  (m/s),  qo  includes  other  sources/sinks  (m/s),  and  kox  and  koy  are  the 

diffusion  coefficients  (m2/s)  determined  by  kax  =  /?5/y(V2Ql/2  j  and  koy  =  h5B /  (n2yQ 1/2 )  with 

1/3 

nx  and  n  being  the  Manning  roughness  coefficients  (s/m  )  in  x  and  y  directions 
(anisotropic  surface  roughness),  respectively,  and  Q  is  the  friction  or  energy  slope  operator 
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potential  gradient,  a  small  number  (e.g.,  1 0  to  10'  )  is  used  for  Q  to  avoid  a  division  by 
zero  in  the  expressions  of  diffusion  coefficients. 
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On  the  upland  area,  raindrops  first  generate  interrill  (sheet)  erosion,  and  then  with 
increasing  flux  overland  flow  may  cause  rill  erosion.  The  soil  erosion  and  transport  by 
overland  flow  can  be  described  by  the  2-D  depth-averaged  NEST  equation  as  described  by 
Equation  (27).  The  bed  change  rate  due  to  size  class  k  is  determined  using  Equation  (28), 
where  Ctk  is  the  depth-averaged  concentration  of  the  k size  class  of  total-load  sediment, 
Dik  is  the  interrill  erosion  rate,  and  D/k  is  the  rill  erosion  rate.  The  bed  material  sorting  is 
simulated  using  the  multiple -layer  approach  presented  in  the  previous  chapter. 
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Interrill  erosion  depends  on  soil  properties,  surface  slope,  vegetation,  land  use,  rainfall 
intensity,  and  hydraulic  factors  of  runoff.  Soil  detachment  by  raindrop  impact  and  splash  is 
the  main  source  of  interrill  erosion.  The  interrill  erosion  rate  is  determined  with  the  Xiang 
(2002)  formula,  which  was  derived  using  the  grey  relation  and  regression  analysis  based  on 
experimental  data  and  is  illustrated  by  Equation  (29),  where  dk  is  the  diameter  of  k  size  soil 
particles,  S0  is  the  bed  slope,  pbk  is  the  bed  material  gradation  in  the  surface  layer,  and  Rc  is 
the  saturated  sediment  transport  capacity  of  interrill.  Low’s  (1989)  formula,  as  given  by 
Equation  (30),  is  used  to  determine  the  saturated  sediment  transport  capacity  of  interrill, 
where  0  =  r  /  [(/?.  - p)gdk  J ,  @c=Tc/[(ps-p)gdk],  r  is  the  flow  shear  stress,  r .  is  the 
critical  shear  stress,  and  Sf  is  the  energy  slope. 
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Sediment  transport  in  rills  is  similar  to  that  in  open  channels.  Based  on  the  non¬ 
equilibrium  concept,  the  rill  erosion  and/or  deposition  by  overland  flow  are  estimated  as 
given  by  Equation  (31).  In  the  equation,  qs,k  is  the  sediment  transport  capacity  of  rill  flow 


calculated  by  Yalin’s  (1963)  formula,  qtk  is  the  actual  sediment  transport  rate,  and  L,  is  the 
adaptation  length  related  to  the  adaptation  coefficient  a,  by  Lt  =  Uhj[atcos0)  . 


Dfk  (9W-  Qtk) 


(31) 


Model  Test  Examples 

The  model  has  been  verified  and  validated  using  analytical  solutions,  laboratory  data  and 
field  data  (He  et  al.,  2009).  One  case  is  the  field  erosion  experiment  of  Barfield  et  al.  (1983) 
on  a  plot  of  4.6  m  x  22.1  m  with  a  slope  of  0.09.  The  bed  material  in  the  experiment  was 
tilled  and  wet  topsoil  with  a  median  grain  size  d50  =  0.06  mm.  The  rainfall  intensity  used  in 
the  simulation  was  6 1  mm/h,  lasting  for  30  min.  The  flow  and  sediment  discharges  measured 
at  the  downstream  end  of  the  plot  were  compared  with  the  model  results.  The  total  simulation 
time  was  40  min  and  the  time  step  was  3  sec.  The  computational  mesh  was  20x10  in  length 
and  width.  A  constant  infiltration  rate  of  4.5  mm/h  was  used  in  the  simulation  due  to  the  lack 
of  soil  properties.  As  shown  in  Figure  14,  one  can  see  that  the  model  provided  a  good 
prediction  for  both  flow  and  sediment  discharges. 
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Figure  14.  Predicted  hydrograph  and  sedimentograph  with  observed  data  of  Barfield  et  al. 

(1983) 


The  other  test  example  was  simulation  of  rainfall-runoff,  soil  erosion  and  transport,  and 
contaminant  transport  in  the  Deep  Hollow  Lake  watershed,  located  in  the  northern 
Mississippi,  U.S.A.  The  watershed  is  about  202  ha,  covered  with  cotton  and  soybean  fields 
and  wood  area.  The  oxbow  lake  has  a  length  of  ~1  km  and  a  width  of  -100  m.  Upland 
elevations  within  the  watershed  range  from  35  to  39  m  above  the  sea  level.  Water  depth  in  the 
lake  ranges  from  0.5  to  2.6  m,  with  the  greatest  depth  in  the  middle.  Two  gauging  stations, 
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DH1  and  DH2,  were  installed  at  two  of  the  inlets  of  the  lake  by  the  U.S.  Geological  Survey. 
Within  the  watershed,  the  surficial  clay  layer  varies  from  2.4  to  4.6  m  while  the  underlying 
silty-sand  layer  is  approximately  33.5  m  thick. 

The  model  was  applied  to  simulate  the  storm  event  of  February  10,  1998.  The 
computational  mesh  size  was  380x252x50  in  length,  width,  and  depth.  The  horizontal  grid 
spacing  (Ax,  Av)  was  5  m,  while  the  vertical  grid  spacing  (A z.)  was  0.2  m.  The  time  step  was 
1.5  sec.  The  soil  mixture  was  divided  into  three  size  classes  of  clay  (0.0017  mm),  silt  (0.027 
mm),  and  sand  (0.805  mm).  The  Manning’s  roughness  coefficient  n  had  values  of  0.1,  0.05, 
0.04,  0.08,  and  0.01  for  forest,  cotton,  soybean,  pasture,  and  water,  respectively.  The 
adaptation  coefficient  at  was  calibrated  as  0.5.  The  predicted  average  sediment  concentration 
during  the  storm  at  Station  DH1  was  1686.1  mg/1,  which  was  very  close  to  the  observed  value 
of  1620.0  mg/1.  Figure  15  shows  the  predicted  and  observed  flow  and  sediment  discharges  at 
Station  DH2  during  the  storm  event.  Although  the  sediment  rate  was  slightly  underestimated 
at  the  beginning  of  rainfall,  the  peak  value  and  trend  were  predicted  generally  well. 


Figure  15.  (left)  the  Deep  Hollow  Lake  Watershed  and  (right)  Flow  and  Sediment  Discharges 

at  DH2  on  02/10/98  (He  et  al.,  2009) 


Transport  of  Cohesive/Non-cohesive  Sediment  Mixtures 

Sediments  in  natural  surface  waters  can  be  classified  as  cohesive  (clay  and  fine  silt)  and 
non-cohesive  (coarse  silt,  sand,  gravel  and  cobble).  Non-cohesive  sediments  transport  in 
disperse  particles,  whereas  cohesive  sediments  transport  mainly  in  floes.  Because  the 
magnitude  of  electrostatic  forces  acting  among  these  particles  is  comparable  to  or  larger  than 
that  of  the  gravitational  forces,  the  fine  particles  may  stick  together  forming  aggregates  (or 
floes)  when  they  collide  due  to  Brownian  motions,  turbulent  mixing,  and/or  differential 
settling.  On  the  other  hand,  during  the  transport  process,  larger  floes  may  be  disaggregated 
into  finer  floes  and  single  particles  due  to  high  shear  or  large  eddy  ejection  and  sweeping. 
Therefore,  the  cohesive  sediment  floes  behave  much  differently  from  the  non-cohesive  single 
particles.  Meanwhile,  the  settled  cohesive  deposits  on  the  bed  may  consolidate  due  to  the 
gravity  and  the  overlying  water  pressure. 
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Modeling  of  either  non-cohesive  or  cohesive  sediment  transport  has  been  widely 
investigated  (Wu,  2007).  However,  how  to  simulate  transport  of  non-cohesive  and  cohesive 
sediment  mixtures  has  been  rarely  studied,  but  have  gained  more  and  more  attention  recently 
(Ziegler  and  Nisbet,  1995;  Lin,  2010).  It  has  been  recognized  that  when  the  fraction  of  fine¬ 
grained  sediments  is  larger  than  about  10%,  sediment  mixtures  consisting  of  cohesive  and 
non-cohesive  particles  may  exhibit  cohesive  properties  and  undergo  very  complicated 
transport  processes.  Flocculation  of  the  cohesive  particles  and  consolidation  of  the  cohesive 
bed  may  significantly  affect  the  erosion,  deposition,  and  transport  of  non-cohesive  particles. 
The  interactions  between  cohesive  and  non-cohesive  particles  should  be  taken  into  account  in 
the  simulation  of  the  mixed  sediment  transport.  A  NEST  modeling  framework  of  mixed 
cohesive  and  non-cohesive  sediments  is  described  below. 

Model  Framework 

Sediment  is  defined  to  be  cohesive  if  its  diameter  is  less  than  about  0.01  mm.  The  model 
divides  the  sediment  mixture  into  a  suitable  number  of  size  classes.  Because  the  size  of  floes 
may  change  with  flow  and  sediment  conditions,  the  sediment  size  classes  are  defined  based 
on  the  gradation  of  the  dispersed  sediment  particles  rather  than  the  floe  gradation.  In  addition, 
as  it  is  difficult  to  solve  the  details  of  aggregation  and  disaggregation  processes  using  multiple 
size  classes  for  the  floes  and  a  large  number  of  model  parameters  for  each  size  of  cohesive 
sediment  need  to  be  evaluated,  only  one  size  class  is  usually  used  to  represent  all  cohesive 
particles  in  the  mixture. 

The  transport  of  cohesive  sediment  is  usually  treated  as  suspended  load  and  governed  by 
Equation  (15).  The  sediment  deposition  rate  is  determined  using  Equation  (32),  where  cosf  k  is 

the  settling  velocity,  and  Ctk  is  the  section-averaged  sediment  concentration.  For  non- 
cohesive  sediment,  the  coefficient  ak  is  called  the  adaptation  coefficient,  which  is 
determined  using  the  formula  of  Armanini  and  di  Silvio  (1988)  or  calibrated  by  measurement 
data.  For  cohesive  sediment,  ak  is  called  the  deposition  probability  coefficient  (Krone, 
1962).  It  is  determined  by  the  formula  of  Mehta  and  Partheniades  (1975)  as  given  by 
Equation  (33),  where  rhd  mjn  is  the  critical  bed  shear  stress  below  which  all  sediments  are 

deposited  on  the  bed,  rbd  max  is  the  critical  bed  shear  stress  above  which  all  sediments  remain 
in  suspension  yielding  a  zero  deposition  rate. 
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For  non-cohesive  sediment,  the  settling  velocity  cosf  k  is  determined  using  a  formula  of 

single  particle  settling  velocity,  such  as  Wu  and  Wang’s  (2006)  formula.  For  cohesive 
sediment,  Thorn  (1981)  and  others  established  several  empirical  formulas  for  the  settling 
velocity  of  the  floes.  Generally,  the  settling  velocity  of  the  floes  can  be  determined  with  the 
relation  given  by  Wu  and  Wang  (2004b)  as  described  in  Equation  (34),  where  cof  is  the 
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representative  settling  velocity  of  the  floes,  cod  is  the  settling  velocity  of  single  particles  of 
the  cohesive  sediment,  Kd ,  Ks ,  Km  and  Kt  are  the  correction  factors  accounting  for  the 
influences  of  sediment  size,  sediment  concentration,  salinity,  and  turbulence  intensity, 
respectively. 

—  =  KdKsKsaK,  (34) 

°\< 

Flocculation  intensifies  as  the  sediment  size  decreases.  Following  Migniot  (1968),  Qian 
(1980),  Huang  (1981),  and  Dixit  et  al.  (1982),  the  correction  factor  Kd  is  evaluated  as  given 
by  Equation  (35),  where  dk  is  the  diameter  of  the  dispersed  cohesive  sediment  particles,  dr 
is  a  reference  diameter,  and  nd  is  an  empirical  exponent.  The  reference  diameter  dr  has  a 
range  from  0.011  to  0.022  mm  based  on  measured  data,  and  nd  is  evaluated  between  1.8  and 
2.0.  Equation  (35)  is  only  suitable  in  the  range  of  dk  <  dr .  For  dk  >dr,  Kd  is  set  to  1.0, 
which  means  no  flocculation  occurs  for  coarse  sediments. 


Kd  =  (dr  /  dk  r  (35) 

As  sediment  concentration  increases,  the  floe  settling  velocity  increases  until  it  reaches  a 
maximum  value  and  then  decreases.  In  analogy  to  the  method  of  Thorn  (1981),  Ks  is 
evaluated  by  Equation  (36),  where  C  is  the  sediment  concentration  in  kg/m  ,  n  is  an 
exponent  ranging  from  1  to  2  with  a  mean  value  of  1.3,  r  is  an  exponent  ranging  3  to  5,  k] 

and  k2  are  coefficients,  k  is  equal  to(l  +  k1C“')/(l  —  k2CpJ  ,  and  Cp  is  the  sediment 

concentration  at  which  the  maximum  settling  velocity  reaches  and  varies  from  1.5  to  15 
kg/m  depending  on  the  sediment  and  water  properties. 

fl  +  k.Cn  0  <C<C„ 
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As  salinity  increases,  the  floe  settling  velocity  increases  and  approaches  to  a  maximum 
value.  The  trend  can  be  approximated  by  Equation  (37),  where  nsa  is  an  empirical  exponent, 


Csap  is  the  salinity  at  which  the  influence  of  salinity  tends  to  be  saturated,  and  Csa  mm  is  a 
small  threshold  value  of  salinity  above  which  Equation  (37)  is  valid.  According  to  the 
investigations  of  Owen  (1970),  Huang  (1981),  and  Yue  (1983),  Csap  has  a  value  of  about  30 


ppt.  The  exponent  nsa  is  0.53  based  on  Huang  (1981)  data,  0.75  according  to  Yue  (1983),  and 
0.03  according  to  Peng  (1989).  This  difference  may  be  because  their  experiments  are  in 
different  ranges  of  salinity. 
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As  the  flow  turbulence  increases,  flocculation  intensifies  first  and  then  reduces.  By  using 
the  experiment  of  Haralampides  et  al.  (2003),  Kt  is  evaluated  using  Equation  (38),  where  kn , 


nn ,  and  nt2  are  empirical  coefficients,  and  rp  is  the  threshold  bed  shear  stress  at  which  Kt 


has  the  maximum  value.  Like  many  other  parameters,  vp  depends  on  the  sediment  properties 
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and  flow  conditions  and  is  about  0.17  N/m"  in  the  experiments  of  Haralampides  et  al.  (2003). 
The  parameters  kn ,  nn ,  and  nt2  need  to  be  evaluated  for  different  cases. 

q<t-tp 

K,=\  (38) 

(1  +  ^i)(r6  / TP)  ”'2  r>r„ 

The  sediment  erosion  rate  is  determined  by  Equation  (39),  where  phk  is  the  fraction  of 
the  k'h  size  class  in  the  surface  layer  of  bed  material,  and  Ebk  is  the  potential  erosion  rate  of 
the  k"'  size  class.  For  non-cohesive  sediment,  Ebk  is  determined  using  Equation  (40),  where 
C*k  is  the  potential  transport  capacity  in  terms  of  concentration  of  bed-material  load  and  is 
determined  using  the  formula  of  Wu  et  al.  (2000). 

Ebk  =  PtkKk  (39) 

Ebk  =  a  k^sf.k^tk  (4^) 

Cohesive  sediment  is  eroded  from  the  bed  either  in  particles  (surface  erosion)  or  in  layers 
(mass  erosion).  The  surface  erosion  occurs  when  the  applied  shear  stress  exceeds  a  certain 
critical  shear  stress,  and  mass  erosion  happens  when  the  applied  shear  stress  exceeds  the  bulk 
strength  of  the  sediment.  Partheniades  (1965)  found  that  the  surface  erosion  rate  is  a  linear 
function  of  the  dimensionless  excess  shear  stress,  while  Raudkivi  and  Hutchison  (1974) 
observed  that  it  is  exponentially  proportional  to  the  dimensionless  excess  shear  stress.  Mehta 
(1986)  clarified  that  the  exponential  relation  is  valid  for  partially  consolidated  beds  and  the 
linear  relationship  is  valid  for  fully  consolidated  beds.  For  cohesive  sediment,  E*bk  is 
determined  by  Equation  (41),  where  M  is  a  coefficient  related  to  the  material  properties 
(such  as  mineral  composition,  organic  material,  and  salinity),  n  is  an  empirical  exponent,  and 
Tce  is  the  critical  shear  stress  for  erosion.  The  critical  shear  stress  for  erosion,  rce ,  depends  on 
the  mud  dry  density,  organic  material,  temperature,  pH  value,  the  Sodium  Absorption  Ratio 
(SAR),  etc.  Partheniades  (1965)  found  n  to  be  1.  In  the  study  of  fine  sediment  transport  in 
lakes  and  reservoirs,  Gailani  et  al.  (1991)  found  that  the  power  index  n  ranges  from  2  to  3 
and  the  coefficient  M  is  related  to  deposition  time. 

Y 

Tb  Tce  (41) 

7 ce  y 

The  critical  bed  shear  stress  for  surface  erosion  is  significantly  affected  by  consolidation. 
It  increases  as  the  dry  bed  density  increases  with  time  and  depth.  Thus,  the  cohesive  sediment 
is  more  difficult  to  erode  further  down  the  bed  layer.  Nicholson  and  O’Connor’s  (1986) 
formula,  Equation  (42),  is  used  to  determine  Tce  as  a  function  of  dry  bed  density,  where  pd0 
is  the  initial  dry  bed  density  (kg/m  ),  rceQ  is  the  initial  critical  shear  stress  (N/m")  at  pd0 ,  kr 
is  a  coefficient  with  a  value  of  0.00037,  and  nT  is  an  exponent  given  a  value  of  about  1.5. 

Tce=TceO+K(Pd-PdoT  (42) 

When  non-cohesive  and  cohesive  sediments  coexist  in  the  bed  material,  the  incipient 
motion  of  non-cohesive  sediment  may  be  affected  by  the  cohesion  if  the  fraction  of  cohesive 
sediment  is  larger  than  a  certain  value  (about  10%).  If  the  fraction  of  cohesive  sediment  is 
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smaller  than  this  value,  denoted  as  pcmin ,  the  model  sets  the  critical  bed  shear  stress  for  non- 
cohesive  sediment  the  same  as  that  when  no  cohesive  sediment  exists.  If  the  cohesive 
sediment  is  dominant  in  the  bed,  say,  the  fraction  reaches  to  a  certain  amount  (about  50%) 
denoted  as  prnm  ,  the  sediment  mixture  is  considered  as  fully  cohesive  sediment,  i.e.,  the 
critical  bed  shear  stress  of  non-cohesive  sediment  is  equal  to  that  of  cohesive  sediment.  In  the 
case  of  cohesive  sediment  fraction  between  p  .  and  n  ,  a  linear  relation  between  the 
critical  shear  stress  and  the  cohesive  sediment  fraction  is  assumed.  Thus,  the  critical  shear 
stress  of  non-cohesive  sediment  size  classes  is  determined  using  Equation  (43),  where  vck  n  is 


the  critical  bed  shear  stress  of  the  k"'  size  class  in  the  situation  where  only  non-cohesive 
sediment  exists,  calculated  using  the  method  of  Wu  et  al.  (2000)  that  considers  the  hiding  and 
exposure  effects  of  non-uniform  bed  material. 


=  i  r. 


v  ck,n 


k,n  ce  ^ ck,n  ) 


Pc  ~  Pc 


c  max  r  c  min 


Pc  <  Pern in 
Pc  mill  ^  Pc  £  Pc 
Pc  >  Pc  max 


(43) 


Equation  (43)  considers  the  interactions  between  cohesive  and  non-cohesive  bed 
sediments.  It  should  be  noted  that  this  approach  has  simplifications  and  assumptions.  More 
sophisticated  approach  may  be  developed.  However,  the  key  issue  is  that  there  is  very  little 
data  in  the  literature  regarding  such  interactions.  Laboratory  experiments  and  field 
investigations  are  needed  to  improve  the  understanding  and  modeling  of  this  complex 
phenomenon. 

The  rate  of  change  in  bed  sediment  mass  and  the  bed  material  sorting  are  simulated  using 
the  multiple -layer  approach  presented  in  the  previous  chapter.  In  addition  to  bed  change 
caused  by  sediment  transport,  the  bed  lowers  due  to  consolidation.  The  detailed  consolidation 
process  can  be  described  by  the  theory  of  Gibson  et  al.  (1967).  The  model  based  on  this 
theory  calculates  the  evolution  of  the  void  ratio  of  a  soil  layer  using  a  1-D  unsteady 
differential  equation  in  the  vertical  direction.  This  approach  is  usually  time-consuming. 
Therefore,  a  simpler  approach  is  often  used,  which  uses  empirical  functions  to  determine  the 
evolution  of  dry  density.  The  dry  bed  density  increases  very  rapidly  in  the  first  a  few  days  of 
consolidation,  then  increases  slowly,  and  finally  approaches  to  the  ultimate  mean  bed  density. 
The  model  approximates  the  temporal  variation  of  dry  bed  density  in  the  first  year  of 
consolidation  using  the  modified  Hayter  (1983)  formula  as  given  by  Equation  (44),  where 
pdl  is  the  dry  bed  density  at  one -year  consolidation  time,  t  is  the  consolidation  time,  and  ap 


and  p  are  empirical  coefficients. 


(44) 

Pci 

When  the  consolidation  time  is  larger  than  (approximately)  one  year,  the  consolidation 
process  is  very  slow.  This  long-term  temporal  variation  of  dry  bed  density  is  calculated  using 
the  Lane  and  Koelzer  (1953)  formula  given  by  Equation  (45),  where  (5  is  a  coefficient. 
Equations  (44)  and  (45)  are  applied  to  only  cohesive  sediment  size  classes.  For  the  entire  bed 
material  mixture,  the  dry  density  is  determined  by  weight-averaging  all  size  classes’  values 
using  the  Colby  (1963)  method. 
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Pd=Pdl+^ogt  (45) 

To  more  accurately  account  for  the  non-uniformity  of  bed  material  in  depth,  multiple 
layers  can  be  used  to  represent  the  bed  material.  The  change  of  each  layer  thickness  is 
determined  using  Equation  (46),  where  t)';  is  the  thickness,  and  pdj  is  the  dry  density  of  the 


j'h  layer  of  bed  material.  Discretization  of  Equation  (46)  and  summation  over  all  bed  layers 


leads  to  the  overall  bed  change  due  to  consolidation  as  given  by  Equation  (47),  where  J  is 
the  total  number  of  the  divided  bed  material  layers. 


^(V,M 


4v.=Z(‘s;*,-'5;)=l>; 


7=1 
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n+1 

\Pdj 


-1 


(46) 

(47) 


Model  Test  Example 

The  above  model  framework  has  been  implemented  in  the  CCHE1D  channel  network 
modeling  package  (Lin,  2010).  The  resulting  sediment  transport  model  has  been  validated 
using  several  laboratory  and  field  cases.  One  of  the  field  cases  was  the  Lower  Fox  River  in 
Wisconsin,  U.S.  Its  bed  material  contains  a  relatively  large  percentage  of  cohesive  sediments, 
which  accumulate  toxic  organic  chemicals  and  heavy  metals  as  a  potential  contaminant 
source  to  Green  Bay.  The  study  reach  was  approximately  1 1  km  long,  extended  from  a  dam  at 
the  upstream  to  the  river  mouth  (Figure  16).  The  East  River  which  joins  the  Fox  River 
approximately  2  km  upstream  from  the  river  mouth  was  considered  in  the  simulation,  because 
its  flow  discharge  was  approximately  10%  of  the  main  stream  discharge  (Gailani  et  al.  1991). 
A  total  of  14  cross-sections  were  assigned  to  the  computational  domain,  12  for  the  main 
stream  and  2  for  the  tributary.  The  sediment  transport  for  a  10-year  storm  event  from  May  22 
to  June  20,  1989  was  simulated.  Figure  17  shows  the  flow  discharges  and  sediment 
concentrations  at  the  DePere  Dam  measured  by  the  USGS  and  daily  water  level  measured  by 
the  National  Oceanic  and  Atmospheric  Administration  (NOAA),  which  were  used  as 
boundary  conditions  in  the  simulation.  The  time  step  was  15  min.  The  Manning’s  n 
coefficient  was  0.02.  Three  size  classes  of  0.00316,  0.0316,  and  0.447  mm  in  diameter  were 
used  to  represent  the  fine,  medium,  and  coarse  sediments,  respectively. 

The  bed  material  was  divided  into  three  layers.  The  mixing  layer  (first  layer)  and  the 
second  layer  had  the  same  initial  thickness.  The  top  two  layers  were  assumed  to  constitute  a 
1-year  deposit  and  the  bottom  layer  was  assumed  to  reflect  a  10-year  deposit.  Various  values 
were  used  for  the  mixing  layer  thickness,  and  a  value  0.002  m  was  found  adequate.  The 
model  parameters  are  shown  in  Table  1.  Figure  18  compares  the  measured  and  calculated 
sediment  concentrations  at  the  river  mouth.  The  simulation  shows  that  the  calculated 
sediment  concentration  reached  the  peak  when  the  flow  was  moderate  in  the  early  time.  The 
sediment  concentration  started  to  decrease  even  though  the  flow  continued  to  increase  or  was 
kept  approximately  the  same,  because  the  top  two  layers  were  eroded  away  in  the  early  days 
and  the  bottom  layer  was  more  difficult  to  erode.  The  magnitude  of  the  peak  sediment 
concentration  and  the  time  to  peak  were  well  reproduced  by  the  numerical  model. 
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Figure  16.  Sketch  of  the  lower  Fox  River 
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Figure  17.  Flow  discharge  and  sediment  concentration  at  Depere  Dam  and  water  level  at  the 
mouth  of  the  lower  Fox  River  from  May  22  to  June  20,  1989 


Figure  18.  Sediment  concentrations  at  the  mouth  of  the  lower  Fox  River  (Lin,  2010) 
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Table  1.  Model  parameters  used  in  the  lower  Fox  River  case 


Kd 

dr  =0.012  mm,  nd  =1.8 

Settling  velocity  cosf  k 

Ks 

k}  =50.0,  k2  =0.008,  n  =  1.3,  r=4.0,  Cp  =3.5  kg/m3 

K, 

Tp  =0.17  N/m2,  nn  =0.165,  nl2  =0.16,  kn  =1.5 

Ksa 

1 

Deposition  probability  a 

Tw.mm  =0.0  N/m2,  rMmax  =  vce 

Potential  erosion  rate  E*bk 

A7=5.0xl0 3  m/s,  n—2.5,  rce  by  Equation  (42) 
with  rce0  =  0.01  N/m2,  kr  =0.00037,  and  nT  =1.5 

Dry  bed  density 

Equation  (45)  with  pdl  =  1489.0,  1041.0  and  480  kg/m3 
(j  =0.0,  91.0  and  256.0  for  sand,  silt  and  clay,  respectively 

Conclusion 

The  NEST  model  has  been  extended  and  applied  to  simulate  non-cohesive  sediment 
transport  induced  by  rapidly-varying  transient  flows,  by  coastal  currents  and  waves,  in 
vegetated  water  bodies,  and  by  overland  flow.  In  the  case  of  rapidly-varying  transient  flows, 
such  as  dam/levee  break  flows,  storm  surge  and  tsunami  waves,  the  flows  may  be 
transcritical,  the  sediment  transport  is  strong  and  the  bed  changes  rapidly.  Thus,  the  presented 
model  adopts  the  generalized  shallow  water  flow  equations  that  consider  interactions  between 
flow,  sediment  transport  and  bed  change,  as  well  as  shock-capturing  schemes  to  handle  the 
mixed-regime  flows.  In  the  coastal  context,  the  flow  model  adopts  the  phase-averaged 
shallow  water  flow  equations  with  wave-induced  radiation  stresses  and  is  coupled  with  a 
spectral  wave  transformation  model.  The  accompanying  sediment  transport  model  accounts 
for  sediment  entrainment  and  mixing  (diffusion  and  dispersion)  by  currents  and  waves.  In 
vegetated  channels  and  marshlands,  the  flow  model  can  consider  the  porosity  effect  of 
vegetation  and  takes  into  account  the  vegetation  drag  and  inertia  forces  in  the  momentum 
equations.  The  associated  sediment  transport  model  includes  sediment  transport  capacity 
affected  by  vegetation.  For  soil  erosion  and  transport  in  upland,  a  2-D  diffusion  wave  model 
is  used  to  simulate  the  overland  flow.  The  model  considers  the  rainfall,  infiltration,  and 
evaporation,  and  the  related  sediment  transport  model  takes  into  account  the  rill/interrill 
erosion  due  to  raindrop  splash  and  hydraulic  shear.  Even  though  different  flow  models  are 
used  in  these  cases,  the  sediment  transport  models  have  similar  model  formulations,  which 
are  extensions  of  the  general  NEST  model  described  in  the  previous  chapter.  The  differences 
in  these  models  arise  from  evaluation  of  sediment  entrainment,  adaptation  length,  and 
diffusion  coefficient. 

A  general  NEST  model  framework  has  been  also  developed  for  simulating  transport  of 
cohesive/non-cohesive  sediment  mixtures.  It  takes  into  account  the  effects  of  cohesive 
sediment  flocculation,  bed  consolidation  and  interactions  between  cohesive  and  non-cohesive 
bed  materials.  Flocculation  is  considered  by  determining  the  floe  settling  velocity  as  function 
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of  sediment  size,  concentration,  salinity  and  turbulent  intensity.  The  critical  shear  stress  of 
cohesive  sediment  erosion  is  related  to  the  dry  bed  density,  which  varies  with  time  due  to 
consolidation.  The  critical  shear  stress  of  non-cohesive  sediments  is  corrected  when  the 
cohesive  fraction  exceeds  a  certain  limit,  such  as  10%.  The  mixture  transport  is  described 
using  the  same  transport  equation  with  unified  formulas  for  the  erosion  and  deposition  rates 
of  cohesive  and  non-cohesive  sediment  size  classes.  The  model  framework  has  been 
implemented  and  tested  in  a  1-D  channel  network  model  and  can  be  extended  to  2-D  and  3-D 
cases  under  various  flow  conditions. 

It  is  well  accepted  that  sediment  transport  modeling  is  highly  empirical  and  many  model 
parameters  are  case  dependent.  Testing  using  laboratory  and  field  measurement  data  has 
demonstrated  that  the  most  important  parameters  in  the  NEST  models  are  the  sediment 
transport  capacity  (or  entrainment)  formula  and  the  adaption  length.  Among  these  two 
parameters,  the  sediment  transport  (or  entrainment)  capacity  formula  is  more  important, 
because  it  affects  the  simulated  bed  morphology  changes  both  quantitatively  and 
qualitatively.  The  adaptation  length  usually  affects  the  simulated  bed  changes  only 
quantitatively.  More  laboratory  experiments  and  field  measurements  are  needed  to  improve 
the  models  in  different  applications. 
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